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This article presents an on-line tool and its accompanying software resources for the numerical 
solution of basic radiation transfer out of local thermodynamic equilibrium (LTE). State-of-the-art 
stationary iterative methods such as Accelerated A-Iteration and Gauss-Seidel schemes, using a short 
characteristics-based formal solver are used. We also comment on typical numerical experiments 
associated to the basic non-LTE radiation problem. These resources are intended for the largest 
use and benefit, in support to more classical radiation transfer lectures usually given at the Master 
level. 


I. INTRODUCTION 

The theory of radiation transfer is of paramount impor¬ 
tance for astrophysics. Indeed, except for those objects 
in the Solar System close enough to us for being explored 
in situ, from the collection of lunar samples back in the 
70’s to the spectacular landing of the Philae spacecraft 
and its instruments on-board, on comet Ghuryumov- 
Gerasimenko in November 2014, our knowledge of the 
Universe overwhelmingly comes from the analysis of the 
light we collect from distant objects. 

Eor the emblematic case of stars, photons are gener¬ 
ated in the central parts of the body. They first scatter 
across its internal and still opaque layers, finally escaping 
the star at the bottom of its atmosphere or photosphere. 
These photons will continue to be scattered through 
the most external layers (and any possible circumstel- 
lar structures therein, like for instance solar prominences 
lying in the corona) of the star, before reaching the inter¬ 
stellar medium and travel into space down to our instru¬ 
ments. Einally, stellar light will scatter again through 
the Earth atmosphere, when one considers ground-based 
astronomical observations. 

The issue of how radiation transfers along these media 
of distinct physical nature, in terms of density, tempera¬ 
ture, dynamics, magnetic field etc. thus appears as quite 
obvious. And even though we shall discard hereafter any 
discussion about terrestrial atmospheric effects, radiation 
transfer through stellar atmospheres still is a very dif¬ 
ficult problem of physics. It relies indeed on eomplex 
non-linear light-matter interactions (see e.g., Hubeny & 
Mihalas 20l4l^. 

The equation of radiative transfer is very likely to be 
present in all astrophysics courses, more likely at the 
Master level. Analytical solutions are very few, but they 
can quite easily be taught and fully derived within a 
few lectures of introduction to the radiation transfer the¬ 
ory. Despite the often very crude approximations used 


in these cases, such solutions may still be very useful to 
any astronomer willing to “clutch at straws” while facing 
a problem of interpretation of data (or of numerical re¬ 
sults) involving some more or less complicated radiative 
modelling. 

At first glance, the radiative transfer equation (here¬ 
after RTE) appears as a deceptively simple first order 
ordinary differential equation. This is indeed the case 
when the so-called souree funetion is already known, 
and the process of deriving the radiation field from a 
known source function is refered to as the formal solu¬ 
tion. However, in the more general case, the RTE is 
integro-differential, because the source function depends 
on integral terms involving the radiation field itself. The 
general problem of defining self-consistently the radiation 
field i.e., the specific intensity and its first moments, to¬ 
gether with the detailed excitation and ionization states 
of an atmosphere (and therefore the opaeity, as well as 
the spatial variations of the source function), consider¬ 
ing the highly non-linear light-matter interactions which 
usually take place in astrophysical plasmas is definitely, 
and still, a (v ery) dijfieult problem (see e.g., Rutily & 
Ghevallier 200d^. 

A considerable simplification of the problem is brought 
by the assumption of considering that the astrophysical 
plasma permeatted by radiation is in the physical condi¬ 
tions of the so-called Local Thermodynamic Equilibrium 
(hereafter LTE). In such a case, velocity distributions of 
particles follow a Maxwell-Boltzmann distribution, exci¬ 
tation and ionisation stages of every atom (or molecule) 
constituting the plasma follow Boltzmann (for excitation 
equilibrium) and Saha (for ionisation equilibrium) statis¬ 
tics, and the self-consistent source functions are, accord¬ 
ingly, characterized locally by Planck functions (see e.g., 
Ghandrasekhar 196Cp^. More insight about the assump¬ 
tions underlying LTE - and therefore how departures 
from LTE can take place in a stellar atmosphere - can be 
found in the monographs of Hubeny & Mihalas (2014)P^ 
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and Oxenius (1986]P^. One should also realize that re¬ 
alistic radiative modelling, even assuming LTE, is such 
complex already that it actually requires numerical mod¬ 
elling. But even though LTE may be suitable for several 
astrophysical “objects” (e.g., stellar photospheres), de¬ 
partures from LTE should always be considered as the 
most general situation. These effects have indeed been 
identified and studied as early as in the late 60’s, with the 
advent of “numerical radiation transfer” (see e.g., Cuny 
196lP, Mihalas & Auer 196^^^. However a few analyt¬ 
ical solutions to the non-LTE radiation transfer problem 
may be derived and used for the sake of validating any 
numerical approach to the problem. 

Excluding probabilistic methods such as Monte-Carlo 
(see e.g., Auer 196^, Bernes 197^^ and Whitney 2011^ 
for a recent review), we may consider that solutions of the 
non-LTE radiation transfer equation fall into two main 
classes being either difference equation methods (e.g., 
Mihalas & Auer 1969), or stationar'^ iterative meth¬ 
ods (e.g., Olson, Auer & Buchler 1986 and references 
therein). Hereafter we shall focus on iterative meth¬ 
ods, from the so-called A-iteration (or Picard) method 
to the Jacobi-like approximate or Accelerated A-Iteration 
(also known as ALI) method and, finally, the more recent 
and much less popular still Gauss-Seidel and Successive 
over-relaxation (SOR) schemes (Trujillo Bueno & Eabi- 
ani Bendicho 1995]P^. 

The treatment of the RTE thus appears to be also a 
good introduction to numerical techniques, including the 
use of the moments of a function, in physical sciences. 

We first remind basic equations driving the non-LTE 
(unpolarized) radiation transfer problem in a static and 
ID plane-parallel geometry. We shall also restrict our¬ 
selves to the special case of monochromatic scattering. 
Then we shall be able to derive an analytical solution of 
the NLTE radiation problem. This solution shall there¬ 
fore be used for testing several iterative methods, includ¬ 
ing the very popular ALL Einally, we shall describe a 
new on-line tool designed for educational purposes. It 
is located at http://rttools.irap.omp.eu/, and the 
associated Python software will also be made available 
from us. 


II. FORMAL SOLUTION OF THE RADIATION 
TRANSFER EQUATION 

The derivation of the RTE in a semi-infinite, plane- 
parallel, static, and ID geometry can be found in several 
classic textbooks (Hubeny & Mihalas 20or in the 
e-book of Rutten 200d^. 

The source function is defined as Su = VuIXu i-e., the 
ratio bewteen monochromatic emissivity and the extinc¬ 
tion or absorption coefficient. The optical depth is de¬ 
fined as dTjy = —Xudz where /i = cos(6>) and dz = fids^ 
assuming that z points in the opposite direction of grav¬ 
ity and that the direction cosine /i defines the orientation 
of a ray vs. the z-axis - see Eig. (1). Then RTE comes 



FIG. 1. Geometry of the ID plane-parallel radiative transfer 
problem. A ray emerges from a semi-infinite atmosphere with 
angle 0 vs. the normal to the surface z. The direction cosine is 
usually called p. Along the ray, we use the spatial coordinate 

s. 


into its more familiar form: 


( 1 ) 

dTy 

where ly and Ty also depend on p. This basic RTE is 
the one which appears in most lectures notes and text 
books introducing radiation transfer to astrophysicists. 
We shall also assume that the source function is isotropic 
i.e., angle-independent. 

The so-called “formal solution” is the general solution 
of this equation for a known source function. In such as 
case, the solution can be easily derived as: 


J ri 

( 2 ) 

We shall see hereafter how this formal solution is used 
together with iterative methods we are particularly in¬ 
terested in. Note also that, when ri ^ 0 and r 2 ^ Too 
in Eq. 0, the formal solution is the Laplace transform 
of the source function. 


A. The Eddington approximation 

It is usual and convenient, for analytical radiation 
transfer, to define the three successive (angular) moments 
of the specific intensity, or Eddington moments: 
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where J is the mean intensity, H the Eddington flux 
(which relates to the astrophysical flux usually observed 
for spatially unresolved objects such as most stars other 
than the Sun), and K which is related to the radiation 
pressure (see e.g., Hubeny & Mihalas 20 

In a similar fashion, successive moments of the RTE 
can be easily derived. It leads, respectively for Eddington 
factors Hjy and to 


dH, 

dry 


Ju-S 
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and 


J(t) = A[5(t)] . (9) 

Eor all computations presented hereafter, as well as 
for the tools we distribute, this process is done using 
the Python module formal which computes the formal 
solution using the so-called short characteristics method 
(Olson & Kunasz 1987^, Auer & Paletou 1994^. 


III. AN ANALYTICAL SOLUTION TO A NLTE 
RADIATION PROBLEM 

An analytical solution to the problem of non-LTE ra¬ 
diative transfer can be derived with the following assump¬ 
tions. Eirst, we shall consider the case of monochromatic 
or coherent scattering. We shall also consider a source 
function that contains a thermal emission component eB 
and a coherent isotropic scattering term (1 — 5 ) J, that is 

S = {l-e)J + eB . (10) 


dK, 

dry 


= H,. 


( 7 ) 


Then, considering the state of the radiation field at 
great depth in a stellar atmosphere when it can be safely 
considered also as isotropic, one can establish the so- 
called Eddington approximation, Jj^ = (Rutten 

2002P^, Hubeny & Mihalas 20l4I^. 

Although valid only in the above-mentioned condi¬ 
tions, it is common in analytical radiation transfer to 
consider that the Eddington approximation remains valid 
throughout the whole atmosphere of a star, even though 
the anisotropy of the radiation field increases while we 
move towards its most external layers. 


B. A-operator 

It is usual in the field of astrophysical radiation transfer 
to write the formal solution of RTE as 


In the frame of the two-level atom model, 5 is also called 
the collisional destruction probability factor (it may also 
be related to a so-called albedo, with a = 1 — e though). 

Assuming that the Eddington approximation J = 3K 
is valid at all depths in the atmosphere, we get easily 
after forming the second derivative of K that 




= 3{J-S). 


(11) 


Eor an isothermal atmosphere of constant B with r and 
using the expression of the source function introduced at 
Eq. (10), the latter expression turns into 


( 12 ) 

whose solution is such that {J — B) = ae~^'^. 

Einally, using the boudary condition at the surface 
J(0) = \/3H{fi) derived by Krook (1955]P^, one can es¬ 
tablish that 


J.(r)=A[5(T)] (8) 

where the operator A represents the operation of deriv¬ 
ing the specific intensity from known opacity and source 
function spatial distributions. It is also, in other terms, 
an integration operator of the known source function 
weighted by the exponential kernel ed~^\ 

Should we ignore the frequency dependence of the ra¬ 
diation field, also known as the “grey case”, and include 
the angular integration leading from specific intensity to 
the mean intensity (i.e., a physical quantity proportional 
to the more generally observed astrophysical flux), one 
would rather write the formal solution as: 


B 


a = — 


1 + V^’ 


(13) 


which leads to the Eddington solution of the non-LTE 
radiation transfer problem: 


S^mJB = 1 - (1 - . (14) 

It is important to identify two critical values associated 
with this solution. Eirst is the surface value, 5'Edd.(^) 
for r ^ 0, of the source function given by ^/eB. Any 
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numerical solution have to be tested vs. this limit value, 
and with some accuracy. Second is the typical depth at 
which 5'Edd. — B which is often called the thermalisation 
depth. It scales as 1/ ^/s for the monochromatic scattering 
case we consider here. Again this typical length should 
be identified with accuracy from any numerical solution 
to the non-LTE problem. 

We conclude this section emphasizing the fact that 
the Eddington solution we have established (for a semi¬ 
infinite atmosphere) is the true solution of an approxi¬ 
mate radiation transfer problem. The main approxima¬ 
tion used here is on the angular dependence of the specific 
intensity. While the resolution of the full problem would 
imply to use an infinite number of directions, we down¬ 
sized it to a single point angular quadrature (see also e.g., 
the discussion provided in §5. of Chevallier et al. 2003). 


IV. NUMERICAL SOLUTIONS 

We shall consider hereafter monochromatic (or grey) 
radiation transfer, so we can drop any frequency de- 
pendance of the Eddington factors in the remainder of 
this article. Boundary conditions are also assumed to be 
monochromatic. 

The Eddington approximation is also compatible with 
the so-called “two-stream approximation”. In that case, 
we adopt also a simplified angular quadrature using 
fi = (i.e., the Van Vleck angle). Beyond astro¬ 

physics, this approximation is also common for global cir¬ 
culation or weather forecasting models developed in (ter¬ 
restrial) atmospheric sciences (J.-P. Chahoureau^ private 
communication). Note again that a proper comparison 
between numerical solutions and the (analytic) Edding¬ 
ton solution requires the use of a single point angular 
quadrature. 

Our formal solver uses short characteristics (SC) us¬ 
ing monotonic parabolic interpolation, as originally de¬ 
scribed in Auer & Paletou (1994 - see also Olson & Ku- 
nasz 1987^, Kunasz & Auer 198^^, and Paletou & Leger 
2007 !^.® Hereafter we only remind the essential princi¬ 
ples of SC, and encourage the reader to consult the exist¬ 
ing scientific literature, for details. The numerical imple¬ 
mentation of short-characteristics rely on the following 
principles. Short characteristics means that the formal 
solution across the whole atmosphere will be carried-out 
depth after depths from one boundary surface to the other 
one, and back-and-forth i.e., for fi negative first then for 
/i positive (note that this order is indifferent, but the sep¬ 
aration between positive and negative direction cosines is 
very important, and shall prove very useful for the case 
of Gauss-Seidel iterations). In order to perform at each 
spatial depth the formal solution expressed by Eq. (2), we 
shall first assume that the source function is quadratie in 
the optical depth. This assumption allows to derive an 
analytical expression of the integral term on the source 
function spatial distribution. Then it can easily be shown 
that, for a current position k the integral in Eq. (2) will 


only involve quantities known at this very position, and 
at the “upwind” or (/c — 1), and “downwind” or (/c + 1) 
positions. Einally, at each depth the current value of 
the specific intensity (for a given direction cosine) will be 
given by 

dk = Ik-1^ : (13) 

where the T’s are analytical functions of the optical 
depths Atu and Ar^ i.e., between the local (0), and the 
upwind (u) and downwind (d) spatial positions. 

Eor all cases discussed hereafter, we use boundary 
conditions more usually used for “semi-infinite” atmo¬ 
spheres. There is no radiation falling (/i < 0) onto the 
top surface of the atmosphere, and the bottom and up¬ 
ward (/i > 0) boundary condition is that the specific in¬ 
tensity equals the Planck function (set to one hereafter). 

A. A— iteration 

The so-called A-iteration (li) is equivalent to a Pi¬ 
card iterative scheme (itself being a fixed-point iterative 
scheme for ODEs). 

Let be the spatial distribution of the source function 
across the atmosphere from (the initial value or) the pre¬ 
vious iteration step. A-iterating consists in successively 
computing 

J = A[5t], (16) 

then = [l — e)J^eB and so on, until convergence. 

Unfortunately, this poor numerical scheme is still 
in use although it is well-known that it is “pseudo- 
convergent” (see e.g., Hubeny & Mihalas 2014^^, Olson et 
al. 1984^21). This is well demonstated by Eigs. (2) where 
we displayed (a) the successive iterates of the Ll scheme 
together with the target analytical solution of Eddington 
and, (b) the respective relative correction, from an iter¬ 
ation to another. Re and the so-called “true error” Tg, 
which is the relative error vs. the analytical Eddington 
solution. 

Eollowing the definitions found in the original papers of 
Auer et al. (1994]Pand Trujillo Bueno & Eabiani Bendi- 
cho (1995/^, these two latter quantities are respectively 

Re = max(| 5 '^^) - |/5'(n)) (17) 

Te = max(|5'^’^^ - 5'Edd. l/5'Edd.) • (18) 

Both are indeed useful to demonstrate the pseudo- 
convergent nature and the failure of Ll. Indeed on Eig. 
(2b) once can notice the constant drop of the relative cor¬ 
rection Re giving the misleading impression that pushing 
the iteration number will finally reach the solution, while 
the true error Tg indicates that the pseudo-solution will 
remain far away from the reference solution of Eddington. 
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FIG. 2. A-iteration for a semi-infinite slab of total optical 
thickness 10^ with 5 points per decade, and an atmosphere 
such diS £ — 10“^. The top hgure shows the successive runs 
of the source function, and the pseudo-convergent nature of 
the numerical scheme. The analytical solution is represented 
by the dotted curve. The bottom hgure displays (i) the maxi¬ 
mum relative correction, from an iteration to another (dashed 
line), and (ii) the true error i.e., the relative error vs. the an¬ 
alytical solution (note that we always have Te > Re after a 
few iterations). 

Note also that, in all hgures > Re after the very 
hrst iterative steps, if not at the onset of the iterative 
process. 


B. ALI: Accelerated/Approximate A— iteration 

The ALI method is basically an operator splitting 
method. Let us write therefore: 


A = A*+M. (19) 

At this point several choices for an approximate opera¬ 
tor A* are possible. However, in our study we shall only 
consider the most efficient version of ALI which is just a 
Jacobi-type method. In such a case. A* should be the 
exact diagonal of the full operator A. The study of ref¬ 
erence concerning this very method is the seminal article 


FIG. 3. Same as Fig. (2) but for Accelerated A-iteration. 
Within the same number of iterations (150) the analytical 
reference solution is reached, unlike with Li. The true error 
reaches a constant level around 0.01, also indicating that there 
is no need to iterate further than ~ 140 iterations, while the 
relative correction continues to decrease. 

of Olson, Auer & Buchler (1986)P3. In practice, the di¬ 
agonal operator is very easily determined, as described 
in Auer & Paletou (1994)1^ 

Now let us write the succession of iterates of the source 
function as S = ^ 6S^ where means the source 

function known at the current iterative step. Now using 
the definition of the A-operator J = A[5'] we can write the 
expression of the source function correction explicitely: 

(55 = [1 - (1 - {(1 - e)A[S^] + eB - 5^} . (20) 

When A* is the diagonal of the full operator, at each 
depth in the atmosphere the increment of source function 
is just obtained from a scalar divide, which make the ALI 
method both accurate and fast. 

Unlike the Ll iteration, it can be rigorously shown that 
the ALI iterative scheme definitely converges to the solu¬ 
tion of the problem, as demonstrated in Fig. (3a). How¬ 
ever, the accuracy of the numerical solution depends on 
the spatial sampling of the atmosphere. It can be mea¬ 
sured by the true error, Tg, which reaches a plateau at 
^ 10“^, as can be seen in Fig. (3b). 
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FIG. 4. Illustration in support of the basic principle of the 
implementation of Gauss-Seidel iterations within the short 
characteristics method. 



FIG. 5. True error and relative correction (dashed line) for 
the same parameters than in previous hgures, but using the 
GS iterative scheme. The truncation error seen by the plateau 
at Te ~ 10is reached in less than 30 iterations now. 


C. Gauss-Seidel and SOR iterations 


Experimenting both Gauss-Seidel (gs) and Successive 
over-relaxation (sor) are logical steps after having ex¬ 
perienced the Jacobi-type ALI methods. Although pub¬ 
lished twenty years ago now, by Trujillo Bueno & Fabiani 
Bendicho (1995)P^, it is not yet of common practice, un¬ 
like ALL 

The GS iterative method consists essentially in updat¬ 
ing the current source function value once the full angu¬ 
lar integration of the specific intensity can be performed 
- because all the necessary quantities yet are available, 
and before the formal solver will be moving to the next 
depth point. This is made relatively easy within the short 
characteristics methods which separates sweping the at¬ 
mosphere first for /i < 0 and second for the remainder 
/i > 0 directions (or vice versa, the important point being 
an explicit distinction between positive and negative di¬ 
rection cosines). This is sketched in Fig. (4) which should 
be read from left to right. Assume that all depths have 
been covered along SC’s of /r < 0, so that we know all spe¬ 
cific intensities and mean intensities for these direction 
cosines, up to the bottom boundary layer N. The next 
task is to complete the angular integration for the up¬ 
ward direction cosine(s). Starting at the bottom bound¬ 
ary, the “upwind” specific intensity is known since it is 
provided by the external boundary condition. Therefore, 
at layer N we can easily compute J^. This knowledge 
makes it possible to update the local source function be¬ 
fore switching to the next inner layer at {N — 1). This is 
the main “trick” for doing Gauss-Seidel iterations within 
the SC method. It requires however modifications of the 
more traditional formal solver used for ALI (provided by 
us as the formalGS module). Indeed, when moving to 
the next layer for /i > 0 at (A^ — 1), we shall advance the 
specific intensity according to Eq. (15) but now using a 


mixture of the just updated and of i.e., of 

the not yet modified values of the source function at the 
local and at the downwind position {N — 2). The process 


is then repeated up to the top boundary surface of the 
atmosphere. 

The numerical implementation of the method was de¬ 
scribed in every detail in the original article of Tru¬ 
jillo Bueno & Eabiani Bendicho (1995)P^. We therefore 
strongly encourage the reader to study this article with 
great care. 

Eigure (5) shows the significative gain on the conver¬ 
gence rate provided by GS. Eor ID problems, it is by far 
superior to the small additional computations induced 
by the indispensable modifications of the classical short 
characteristics formal solver. 

Einally the SOR scheme is built on the same strategy 
adopted for GS although the new source function incre¬ 
ment with uj chosen between 1 and 2. 

It can be shown that the optimal scheme is obtained for 
uj ^ 1.5 (see Trujillo Bueno & Eabiani Bendicho 199^^ 
More insight about the SOR method can be found in 
Young (1971).!^ 


V. NLTE RADIATION TRANSFER ON-LINE 

We have implemented a dedicated web-service which 
allows on-line numerical experiments with the nu¬ 
merical methods we just presented. It is located 
at http: //rttools . irap. omp. eu/ , and maintained by 
OMP-IRAP (Toulouse, France). In addition, a Git repos¬ 
itory is about to be installed, in order to distribute the 
original Python modules we developed and used for the 
web-service. 


A. Description 

Several “buttons” may be independently activated. 
The first choice is to be made among methods i.e., A- 
iteration (LI), Accelerated A-iteration (ALI), and Gauss- 
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Seidel or SOR (GS). Concerning the latter choice, the 
distinction between both schemes will be controled by 
the uj (omega) parameter. Once the method have be se¬ 
lected, the user will have to provide a value for 5 (eps) 
with format 10“^ where p is an integer, the total optical 
thickness of the atmosphere Tmax (taumax), using a for¬ 
mat 10^, and the number of points per decade (npdec) 
used for setting the spatial grid for the computations. Fi¬ 
nally, the number of iterations (niter) to be performed 
is required. 

The output consists in two graphics. The first one pro¬ 
vides the history of the source function (initialized with 
S = B) for the number of iterations required, with the 
corresponding Eddington solution vs. optical depth. The 
second plot displays both true error and relative correc¬ 
tion (from an iteration to another) vs. the number of 
iterations. 


B. Experiments 

As a preamble to this part, we remind here that com¬ 
parison with the analytical solution of Eddington make 
only sense for the case of “effectively thick” slabs, that is 
such that Tmax ^ f^LC casc of monochromatic 

scattering we only consider hereafter. 

A first and indispensable experiment is to realize the 
failure of the A-iteration . The “pseudo-convergence” of 
this method can be noticed by the rapid and continuous 
fall of the relative correction while the “solution” at 
which the process is converging may remain very far from 
the Eddington solution, as indicated by the “true” error 
Tg. Eurther experiments pushing the number of itera¬ 
tions for the same set-up should be done and analysed. 

Using the very same set of parameters, the next step 
is to experiment the benefit of the diagonal operator ALI 
method. The latter should quickly reach the Eddington 
solution with good accuracy, unlike Ll. Other experi¬ 
ments will allow to check the so-called y^-law for the 
surface value of 5'/5, as well as to identify the “thermal- 
isation depth” 1 / ^/e. 

However, even ALI may work at limited accuracy. A 
relevant experiment is to iterate the method with a given 
set of parameters up to reaching a plateau for Tg. The 
latter value gives an indication of the truneation error of 
the method, due to the spatial discretization of the nu¬ 
merical method. All other parameters remaining equals, 
one should then experiment the effects of only modify¬ 
ing the sampling of the slab, by changing the number of 
points per decade parameter. Changes in the limit value 
of Tg, together with the rate of convergence should be in¬ 
vestigated. Note that a detailed study on the accuracy of 
the ALI method was published by Chevallier et al. (2003). 


Einally, we propose to go beyond the ALI- Jacobi 
method with the Gauss-Seidel (gs) and Successive Over- 
Relaxation (sor). Its implementation requires several 
touchy modifications in the original short characteristics 
formal solver which require special attention. We found 
them pretty well documented in the original article, and 
the interested user will get to it by a careful inspection 
of the source code that we also deliver with our web- 
service. Gauss-Seidel and SOR differ only by the choice 
of the relaxation parameter 00 . It should be set to cj = 1 
for performing GS iterations, although it should be picked 
between 1 and 2 for experimenting SOR iterations. It is a 
good exercice to test various values of cj, seeking for an 
optimal scheme. 

A more insider study, requiring to use directly the 
Python code we make available, would be to test the so- 
called smoothing capability of GS/SOR methods which 
plays a crucial role in multi-grid methods (see e.g., Auer 
et al. 199#1) 

VI. CONCLUSION 

We made available a tool very suitable to any astro¬ 
physics Master programme, or for any astronomer or 
physicist willing to start an initiation to state-of-the-art 
numerical radiation transfer. 

It is quite straightforward to upgrade the simple angu¬ 
lar quadrature we used for that study. It would be also 
possible to implement, from our Python formal solver, 
the computation of a more realistic scattering integral, by 
integrating an explicitely z/-dependent mean intensity Jj^ 
weighted by an a priori known absorption profile (Gaus¬ 
sian/Doppler or Voigt). 

Possible evolutions may be, to develop a specific for¬ 
mal solver for the ID spherical problem (see e.g., Auer 
198#1), or to propose a simple multi-level atom version 
following the so-called Multilevel-ALI method originally 
developped by Rybicky & Hummer (1991, see also Pale- 
ton & Leger 2007). Gomparisons with (non-stationary) 
conjugate gradient type method could also be set-up (see 
e.g., Paletou & Anterrieu 200^^. 
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